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Executive  Summary 

Pushing  sensors  and  algorithms  to  the  limit  to  minimize  the  chance  of  overlooking  unexploded 
ordnance  (UXO)  increases  the  chance  that  noise  will  be  misidentified  as  signal  and  money 
wasted  excavating  scrap  metal  or  chucks  of  magnetic  rock.  Approaches  to  improve  the  signal  to 
noise  ratio  typically  follow  one  of  three  tracks:  (1)  development  of  new  sensors  that  are  either 
more  sensitive,  or  less  noisy;  (2)  fusion  of  data  from  multiple  sensors  so  that  the  chances  of 
confounding  all  of  the  det6ectors  simultaneously  is  reduced;  or  (3)  development  of  computer 
algorithms  to  extract  the  signal(s)  from  the  noise.  In  the  case  of  wide-area  surveys  that  use 
methods  such  as  helicopter-borne  magnetometry,  the  battle  against  noise  has  focused  primarily 
on  instrument  noise  (e.g.,  thermal  noise),  platform  noise  (e.g.,  magnetic  noise  created  by  the 
helicopter),  or  interference  created  by  multiple  UXO  targets  in  cluttered  settings.  What  has  been 
largely  ignored  is  the  "noise"  created  by  the  rock  and  soil  that  surrounds  the  buried  ordnance. 

The  goal  of  this  pilot  project  was  to  investigate  a  new  approach  to  the  characterization  and 
simulation  of  geologic  noise  using  multifractal  analysis  that  captures  the  scale-dependent 
variability  arising  from  geologic  heterogeneity  in  different  environments.  By  combining  geologic 
noise  simulations  with  models  for  the  geophysical  signatures  of  UXO,  the  researchers  aim  to 
create  synthetic  datasets  that  can  be  used  both  to  test  and  improve  UXO  discrimination 
algorithms  developed  by  other  researchers,  and  to  develop  more  reliable  estimates  of  the  ratio  of 
false  positives  to  false  negatives  for  a  particular  geologic  environment. 

In  this  report  we  discuss  the  multifractal  methodology  and  its  application  to  three  data  sets: 

Isleta  Pueblo,  NM;  Fort  Ord,  CA;  and  the  Sierra  Army  Depot,  CA,  The  SI  site  at  the  Isleta 
Pueblo,  which  is  underlain  by  volcanic  rock,  typifies  a  UXO  site  with  high  magnetic  noise.  Fort 
Ord  is  underlain  by  sedimentary  rock  that  is  largely  non-magnetic,  but  apparently  magnetic 
sediments  have  washed  into  topographic  drainages,  creating  moderate  levels  of  geologic  noise 
distributed  anisotropically  across  the  site.  The  Sierra  Army  Depot  rests  on  a  thick  sediment-filled 
graben  that  is  low  in  magnetic  minerals  so  represents  a  geologic  environment  with  magnetic 
noise  levels  as  low  as  is  found  anywhere.  Thus,  the  three  data  sets  encompass  a  wide  range  of 
geologic  conditions  likely  to  be  encountered  at  UXO  sites. 

Tests  showed  the  data  to  be  multifractal,  and  our  simulation  results  demonstrate  that  the 
multifractal  methodology  provides  a  versatile  tool  for  researchers  to  experiment  with  new 
detection  and  discrimination  algorithms,  and  could  potentially  be  used  for  QA  assessment  at 
UXO  remediation  sites.  Future  work  might  focus  on  improved  methods  for  characterizing  and 
modeling  anisotropy,  on  incorporating  remnant  magnetization,  and  on  joint  multifractal  modeling 
of  other  geophysical  properties  such  as  electrical  conductivity  and  dielectric  permeability. 
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Introduction 


The  Problem  of  UXO  Detection  and  Discrimination 

Scientific  measurement,  whether  looking  at  spectra  from  distant  galaxies  or  searching  for 
irregularities  in  the  EKG  of  a  heart  patient,  invariably  revolves  around  a  struggle  to  extract  signal 
from  noise,  giving  rise  to  the  twin  detection  problems  of  false  positives  and  false  negatives. 

These  two  errors  are  complementary:  raise  instrument  detection  thresholds  to  avoid  false 
positives,  you  increase  the  chance  of  missing  the  signal,  a  false  negative.  Conversely,  lower  the 
detection  threshold  to  avoid  false  negatives,  you  increase  the  risk  of  misidentifying  noise  as 
signal,  a  false  positive. 

The  latter  problem  of  false  positives  is  particularly  acute  when  seeking  to  map  unexploded 
ordnance  (UXO).  To  minimize  the  danger  of  leaving  ordnance  behind  virtually  every  anomaly 
detected  that  could  conceivably  be  a  munition  item  must  be  excavated.  Consequently,  millions 
of  dollars  are  spent  each  year  digging  up  scrap  metal,  and  in  many  cases,  naturally  magnetic  soil 
and  rock.  Butler  (2003)  describes  a  typical  scenario  of  a  $30  million  remediation  program 
where  76%  of  the  budget  is  expended  on  non-UXO  removal  (Figure  1).  In  extreme  cases,  such 
as  Kaho'olawe,  Hawaii,  of  49,521  anomalies  excavated,  only  3%  were  UXO  (Putnam,  2001). 


Example  Removal  Action  Project: 
Cost  Percentages 

ASSUME 

-  Heavily  Contaminated  5000  Acre  Site 

-  Mortar  and  Artillery 


Non-UXO 

Removal 

76% 


UXO 

Removal  9% 

Vegetation 
Removal  4% 


Total  Cost: 
$30  million 


a.  Typical  ordnance  items  b.  Cost  distribution  fir  UXO  cleanup 

Figure  1:  Clean-up  costs  for  a  hypothetical  UXO  removal  scenario  typical  of  remediation  efforts 
in  the  U.S.  (after  Butler,  2003) 
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While  there  has  been  continued  advancement  in  UXO  detection  (e.g.,  Davis  et  al.,  2005),  in 
recent  years,  research  has  focused  many  on  developing  methods  of  discriminating  UXO  from 
non-UXO  items  before  excavation.  Billings  et  al.  (2002),  for  example,  developed  an  algorithm 
that  automatically  fits  a  dipole-model  to  each  magnetic  anomaly  and  compares  the  extracted 
model  parameters  with  a  library  of  UXO  items.  They  found  that  remnant  magnetization  was  the 
most  effective  predictor  of  UXO,  but  that  some  noise  anomalies  remained  indistinguishable  from 
UXO.  Gamey  (2005)  discusses  the  use  of  3D  magnetic  data  and  “gradient  strings”  to  improve 
magnetic  discrimination.  Sanchez  et  al.  (2006)  discuss  importance  of  going  beyond  a  simple 
dipole  model,  and  suggest  discrimination  can  be  improved  by  the  use  of  higher-order  magnetic 
moments  to  obtain  information  about  the  shape  of  the  anomaly  source.  A  common  thread 
running  through  these  research  efforts  detection  is  the  difficulty  of  UXO  detection  and 
discrimination  in  the  presence  of  geologic  noise.  Butler  (2001)  summarizes:  “Although  spatial 
wavelength  filtering  might  be  used  for  UXO  detection  in  the  presence  of  localized 
geologic  anomalies,  UXO  detection  with  magnetometry  in  a  highly  cluttered  environment  will  be 
extremely  difficult  if  not  impossible.”  Considerable  effort  has  been  expended  characterizing  the 
magnetic  signatures  of  an  enormous  variety  of  UXO.  We  believe  that  further  progress  in 
detection  and  discrimination  depends  on  developing  a  similar  in-depth  understanding  of  geologic 
noise. 

Models  for  Geologic  Noise 

To  be  effective,  both  detection  and  discrimination  algorithms  require  a  mathematical  model  of 
the  background  magnetic  noise  arising  from  heterogeneities  in  the  magnetic  susceptibility  of  the 
underlying  soil  and  bedrock.  Traditionally,  researchers  have  assumed  uncorrelated,  uniform 
random  or  Gaussian  random  noise,  but  this  can  lead  to  highly  optimistic  expectations  because 
natural  geologic  noise  is  more  complex.  Searching  for  a  signal  embedded  in  Gaussian  noise  is 
like  trying  to  listen  to  a  conversation  while  standing  close  to  a  hissing  steam  pipe.  The  noise  is 
annoying,  but  different  enough  from  speech  that  you  can  still  make  out  the  words.  Looking  for 
UXO  in  the  presence  of  geologic  noise  is  much  more  analogous  to  trying  to  follow  a 
conversation  in  at  a  crowded  cocktail  party  where  you  are  surrounded  by  people  trying  to  talk 
over  each  other  -  a  much  more  difficult  detection  problem  because  of  the  similarity  between  the 
signal  and  the  noise. 

For  this  pilot  project  we  focused  on  UXO  detection  and  discrimination  using  magnetic  methods. 
Magnetics  has  a  long  and  successful  history  of  application  to  UXO  mapped  precisely  because  in 
many  environments  the  magnetic  signature  of  UXO,  which  are  largely  made  of  steel,  greatly 
exceeds  the  magnetic  background.  This  is  not  true  of  all  sites,  particularly  those  underlain  by 
rock  and  soils  of  volcanic  origin,  for  example  sites  in  Hawaii  and  New  Mexico. 

Researchers  surveying  a  UXO  site  in  New  Mexico  reported  that  over  39%  of  the  magnetic 
anomalies  excavated  were  not  created  by  UXO,  concluding  that  “localized  zones  of  rock  or  soil 
with  high  magnetic  susceptibility  (‘hot  rock/dirt’)  may  have  contributed  to  the  high  rate  of 
‘no  finds'.  ”  (ORNL,  2004). 

In  a  recent  SERDP  project,  Li  et  al.  (2006)  went  beyond  the  usual  assumption  of  uncorrelated 
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random  noise.  They  used  a  simple  model  for  geologic  noise  that  is  Gaussian,  but  with  a  power 
spectrum  that  decays  with  wavenumber  to  create  a  spatial  correlation,  and  used  random  changes 
in  phase  to  generate  different  noise  realizations.  They  then  used  optimal  Wiener  filtering  to  try 
to  separate  UXO  from  background.  This  approach  to  noise  simulation  is  certainly  more  versatile 
than  random  Gaussian,  but  does  not  capture  either  the  scaling  behavior  that  is  important  for 
wide-area  surveys,  or  an  important  characteristic  of  natural  variability  -  zones  of  low  and  high 
noise,  or  intermittency. 

Experience  has  demonstrated  that  many  geologic  processes  are  scale  independent,  meaning  that 
geologic  properties  and  patterns  appear  similar  on  many  scales.  For  this  reason  when  geologists 
photograph  a  rock  outcrop  they  habitually  include  a  rock  hammer,  lens  cap,  or  some  other 
familiar  object  for  scale  (Figure  2),  otherwise  it  is  often  impossible  to  distinguish  irregularities 
on  the  surface  of  a  rock  specimen  from  a  surface  as  large  as  the  lunar  landscape.  The  scale 
independence  of  geologic  properties  has  important  implications  for  modeling  geologic  noise  for 
UXO  detection  and  discrimination  in  wide  area  surveys.  The  implication  is  that  at  larger  scales, 
larger  natural  geologic  background  anomalies  will  be  more  common.  The  geologic  noise 
information  determined  from  data  collected  on  test  site  a  few  hundred  meters  on  a  side  will  not 
reflect  the  amplitude  variability  likely  to  be  encountered  when  surveying  tens  or  hundreds  of 
square  kilometers. 


Figure  2:  Geologic  patterns  are  often  scale  independent.  Here  a  false  scale  on  the  inset 
photo  ( upper  left )  makes  the  side  of  an  island  appear  to  be  a  small  rock  outcrop. 
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For  wide  area  assessments, where  ultra-low  altitude  magnetic  data  are  collected  over  an  entire 
military  base  or  target  range,  it  is  particularly  important  to  use  a  model  for  geologic  noise  that  are 
robust  at  any  scale,  from  sub-meter  to  tens  of  kilometers.  For  this  reason,  we  chose  to  develop 
geologic  noise  models  based  on  multiffactals. 

Fractals  and  Multifractals 

Fractals  can  be  though  of  as  geometric  entities  or  shapes  that  are  created  from  a  simple 
mathematical  generator  that  repeats  on  every  scale.  A  simple  example  is  the  Koch  snowflake. 
Starting  with  an  equilateral  triangle,  we  use  an  algorithm  that  recursively  replaces  the  central 
third  of  each  line  segment  with  two  smaller  segments  of  equal  length  pointing  outward.  Several 
generations  of  snowflake  are  show  in  Figure  3.  When  this  process  is  continued  indefinitely,  the 
result  is  a  shape  has  an  infinite  perimeter  and  no  preferred  scale.  We  can  zoom  in  any  number  of 
times  on  any  portion  of  the  snowflake  perimeter  and  we  will  see  the  same  pattern.  What  makes 
fractals  more  than  mere  mathematical  curiosities  is  that  many  elements  of  nature  show  this  self¬ 
similarity  over  a  large  range  of  scales. 


The  concept  of  a  fractal  can  be  generalized  from  a  simple,  deterministic,  geometric  entity  such  as 
the  Koch  Snowflake,  to  statistical  fractals,  where  is  is  not  the  exact  shape  that  repeats  on  every 
scale,  but  the  statistical  properties  that  are  self-similar  (Figure  4).  In  a  classic  paper  entitled 
“How  long  is  the  coastline  of  Britain,”  Mandelbrot  discusses  how  the  measurement  of  something 
as  simple  as  the  distance  around  an  island  depends  on  the  length  of  the  ruler.  What  Mandelbrot 
found  was  that  there  is  a  direct  statistical  relationship  between  the  length  of  the  ruler  and  the  total 
measured  coastline.  As  the  ruler  becomes  smaller,  we  are  able  to  follow  each  twist  and  turn  of 
the  coast  more  precisely  and  the  resulting  perimeter  grows.  Imagine  using  a  ruler  so  small  that 
you  could  trace  the  boundary  of  each  sand  grain  along  the  water's  edge  -  the  perimeter  of  Great 
Britain  would  be  millions  of  miles!  This  behavior  is  like  that  of  the  Koch  snowflake:  there  is 
increasing  irregularity  as  you  repeatedly  zoom  in  on  the  boundary.  The  self-similarity  is 
statistical  rather  than  geometric,  but  a  fractal  model  captures  many  of  the  features  of  a  coast,  and 
fractal  generators  can  be  used  to  create  extremely  realistic  models  of  coastlines  at  whatever  scale 
is  required.  By  starting  with  different  generators,  an  enormous  variety  of  fractals  can  be  created, 
and  these  fractals  can  be  used  to  generate  extremely  realistic  models  of  leaves,  clouds, 
landscapes,  and  other  natural  objects  (Figure  4).  Thus,  it  is  reasonable  to  suppose  that  a  fractal 
model  might  also  serve  as  the  basis  for  a  scale-independent  generator  natural  magnetic  noise. 

But  first  we  must  make  one  further  generalization,  from  fractal  to  multifractal. 

To  model  noise  requires  statistical  fractals,  but  ordinary  fractals,  which  are  characterized  by  a 
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single  fractal  dimension,  cannot  capture  a  key  feature 
of  natural  processes  -  intermittency.  When  looking  at 
plots  of  natural  magnetic  data  collected  over  a  large 
area  ,  it  is  not  uncommon  to  observe  regions  that  are 
relatively  smooth  punctuated  by  regions  of  higher 
amplitude  magnetic  anomalies.  Furthermore,  the  larger 
the  scale  of  the  data  set,  the  larger  the  chances  of 
encountering  large- amplitude  anomalies.  Multifractals, 
unlike  ordinary  fractals,  can  capture  this  punctuated 
behavior. 


Figure  4.  Fractal  simulation  of  a  cloud. 

A  multifractal  can  be  thought  of  as  a  a  hierarchy  of 
fractal  sets  each  possessing  a  different  fractal 
dimension  that  depends  on  the  field's  threshold  value 
(Figure  5).  This  can  be  made  clear  with  an  example. 

First,  it  is  important  to  realize  that  fractals  are  like 
regular  Euclidean  objects;  when  intersected  with  a 
plane,  the  resulting  intersection  set  has  the  dimension 
of  the  object  minus  one.  Intersecting  a  sphere  (three- 
dimensional)  with  a  plane  results  in  a  disk  (two- 
dimensional).  Imagine  a  fractal  representation  of 
land  surface  elevation  (topography).  If  you  slice 
through  the  surface  at  a  given  elevation  (i.e.,  a 
horizontal  plane),  the  perimeter  of  the  slice  is  also 
fractal,  with  fractal  dimension  equal  to  that  of  the 
surface  minus  one.  If  the  topography  were  viewed 
as  monofractal  of  dimension  D,  then  the  dimension 
of  the  perimeter  is  D-l,  which  is  independent  of  the 
elevation  at  which  you  cut  the  surface.  This  is  not 
realistic,  as  one  expects  the  perimeter  to  be  more 
sparse  as  the  elevation  is  increased.  For  a  multifractal 
surface,  each  perimeter  will  be  fractal  but  with  a 
different  fractal  dimension  (that  decreases  with 
elevation). 

The  multifractal  model  is  viewed  as  a  means  of 
“distributing”  a  quantity  (e.g.,  land  surface  elevation) 
over  space.  It  has  a  strong  physical  basis,  arising 
naturally  from  the  study  of  turbulence  and  FiZure  5-  The  ^me  fractal  cut  at  different 

atmospheric  dynamics,  modeling  the  cascading  of  thresholds.  For  a  multifractal  the  fractal 
energy  from  large  eddies  into  ever  smaller  ones.  In  dimension  of  the  perimeter  (white  line ) 
recent  years,  numerous  papers  have  been  published  */’«■*<>"  threshold,  but  does  not  for 
demonstrating  that  various  natural  processes  ranging  a  monofracta  ■ 

from  cloud  formation,  to  topography,  to  earthquake  periodicity,  can  be  effectively  modeled  using 


7 


multifractals. 


Project  Goals 

The  goal  of  this  project  is  to  investigate  the  use  of  multifractals  to  model  the  magnetic  noise 
encountered  by  wide-area  aeromagnetic  surveys  of  UXO  sites.  More  realistic  models  of 
geologic  noise  can  be  used  to  assist  researchers  developing  detection  and  discrimination 
algorithms  by  creating  an  electronic  proving  ground  of  any  size  that  accurately  captures  the 
variability  encountered  in  nature  for  a  wide  variety  of  geologic  environments.  Multifractal 
models  of  geologic  noise  could  then  be  used  to  assess  the  likelihood  of  false  positives  and  false 
negatives  as  a  function  of  the  size  of  the  area  to  be  surveyed,  the  survey  parameters,  and  the 
geologic  environment.  In  this  pilot  study  we  have  focused  on  the  multifractal  representation  of 
natural  variations  in  magnetic  susceptibility.  The  methods  we  discuss  may  be  generalized  in  the 
future  to  include  other  geologic  properties,  such  as  electrical  conductivity  and  dielectric 
permeability,  which  control  the  noise  levels  for  time-domain  electromagnetics,  ground 
penetrating  radar,  and  other  geophysical  methods  that  are  employed  to  map  UXO. 

Methods 

Universal  Multi  fractal  Modeling 

Fractals  belong  to  the  category  of  scaling  models  that  rely  on  the  fact  that  the  values  of  the  field 
at  different  scales  are  related  through  a  transformation  that  involves  only  the  scale  ratio  (i.e., 
there  is  no  characteristic  length  scale).  Or  mathematically, 

Pr{<PA> Ay)~  A-C(y)  (D 

where  Pr  refers  to  “probability,”  the  tilde  means  equality  within  slowly  varying  constants,  and 
C(y)  is  the  codimension  of  the  field,  a  function  that  is  non-linear,  increasing,  and  concave  with 
the  order  of  the  singularity,  y.  Or  in  plain  English,  the  probability  of  the  multiff  actal  field 
exceeding  the  threshold  kYat  any  given  scale,  X  (meters,  kilometers,  hundreds  of  kilometers),  is 
proportional  to  that  scale  raised  to  the  codimension  C(y),  so  the  larger  the  scale,  the  greater  the 
chances  of  the  field  exceeding  any  given  threshold,  implying  that  you  find  more  large  noise 
anomalies  as  you  increase  the  scale  of  the  survey. 

Equivalently,  we  can  express  the  scaling  properties  in  terms  of  statistical  moments  (e.g.,  mean, 
variance)  instead  of  probabilities: 

<c Pl)~  \Kiq)  (2) 

where  (cpl)  is  the  qth  moment  of  the  field  (with  q  taking  on  non- integer  values),  the  angle 
brackets  represent  ensemble  averages,  and  K(q)  is  related  to  the  codimension  in  a  Legendre 
transform  (see  Boufadel  et  al.,  2000).  Equation  2  states  that  if  we  know  all  of  the  statistical 
moments  of  a  multifractal  field,  this  is  equivalent  to  knowing  its  probability  distribution.  The 
computation  of  moments  is  easier  than  computation  of  the  probability  distribution  from  the  data. 

Scaling  properties  of  a  geophysical  quantity,  G,  may  be  characterized  using  the  structure  function 
(Monin  and  Yaglom,  1975,  Davis  et  al.,  1994;  Boufadel  et  al.,  2000): 
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(3) 


<I^Gj)=  (|G(*)-G(*+s)|,)«j5(,) 

where  s  is  the  spatial  lag,  the  symbol  ~  means  proportionality  within  slowly  varying  quantities, 
such  as  Log(s)  (Meneveau  and  Sreenivasan,  1987),  and  the  angle  brackets  imply  ensemble 
averages.  The  statistical  behavior  of  the  field  is  governed  by  the  structure  function  exponent, 
^(q).  Equation  (3)  generalizes  the  variogram  (q  =  2)  to  include  lower  and  higher  order  moment. 
This  generalization  is  needed  when  the  process  is  multifractal.  Equation  3  states  that  the 
structure  function  (e.g.,  the  variogram)  is  a  power  law  function  of  the  lag  (scale).  This  is  key, 
because  it  tells  us  that  ifc(q)  is  known  at  any  scale  we  can  fully  characterize  G  at  any  other 
scale. 

To  determine  whether  or  not  a  field  is  scaling,  one  computes  <AGs(q)  >  for  a  certain  value  of  q 
(say  q  =  2)  as  a  function  of  the  lag,  s,  and  plots  Log(<AGs(q)  >)  versus  Log(s).  If  the  result  is 
approximately  linear  then  the  quantity  G  is  scaling.  The  slope  of  this  line  yields  q(q)  for  the 
value  of  q  plotted.  To  distinguish  between  monofractal  and  multifractal  behavior,  the  process  is 
repeated  for  many  values  of  q.  One  then  plots  of  c,(q)  as  a  function  of  q.  If  q(q)  varies  linearly 
with  q  (say  up  to  q  <  4)  then  G  is  monofractal.  If  q(q)  is  nonlinear  convex  (downward  facing) 
then  G  is  multifractal.  The  structure  function,  ^(q),  behaves  this  way  because  there  exists  a  one- 
to-one  relationship  between  between  ^(q)  as  a  function  of  q,  and  the  fractal  dimension  as  a 
function  of  threshold  D(T)  (Frisch  and  Parisi,  1985;  Menevearu  and  Sreenivasan,  1987).  Recall 
that  for  monofractals  D(T)  is  independent  of  the  threshold  value  T,  hence  q(q)  increases  linearly 
with  q,  but  such  is  not  the  case  for  multifractals.  Higher  order  moments  (large  values  of  q) 
accentuate  the  role  of  large  values  in  the  field  G  (higher  thresholds),  while  for  smaller  values  of 
q,  c,(q)  depends  primarily  on  the  smaller  values  of  G. 

The  Fourier  power  spectral  density  or  simply  the  spectrum  is  another  tool  to  assess  scaling.  A 
field  G  is  scaling  if  its  Fourier  spectrum  is  power  law  as  function 
of  the  wave  number,  viz: 

EGcck-p  (4) 

where  k  is  the  spatial  wave  number.  The  field  G  is  called  fractal  if  p  is  between  1  and  3;  this  is  to 
allow  D(T)  to  be  less  than  2  if  the  data  are  a  one-dimensional  series.  Such  a  field  is  non¬ 
stationary  and  its  variogram  keeps  increasing  with  the  lag  distance,  s.  For  P  >  3,  the  field  is 
scaling  but  not  fractal.  An  example  would  be  the  distribution  of  wave  heights  at  sea  as  a  function 
of  space,  where  4  <  p  <  5  (Phillips,  1985).  If  p  is  such  that  -1  <  P  <  1,  then  the  field  G  represents 
the  increments  of  a  fractal  field.  It  is  also  stationary  (Davis  et  ah,  1994),  and  the  variogram  of  G 
would  exhibit  a  sill.  We  are  referring  to  the  so-called  “wide-sense  stationarity”,  which  relates  to 
the  behavior  of  the  second  order  moments,  such  as  the  correlation  function  and  the  variogram 
(Yaglom  1987,  p56).  We  believe  that  the  stationarity  of  the  probability  distribution  function  (also 
known  as  strict  stationarity)  is  rarely  achieved  with  geophysical  data.  A  more  recent 
classification  based  on  the  value  of  p  is  the  following:  if  p  <  1  the  field  is  conservative,  for  p  >  1 
(but  still  P  <  3)  the  field  is  nonconservative.  Conservative  fields  have  their  mean  value 
conserved  when  the  scale  of  observation  is  changed. 
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From  an  analysis  point  of  view,  when  P  >  1,  the  structure  function  (SF)  is  a  convenient  tool  for 
characterizing  scaling  properties.  For  fields  that  are  scaling  and  conservative  (i.e.,  -1<  P  <  1),  the 
SF  is  not  a  power  law  (i.e.,  the  second  equality  in  Eq.  3  is  not  valid)  but  the  correlation  function 
is  a  power  law  at  large  lags.  The  correlation  function  is  defined  as: 

R{s)=  ({G{x)  —  G)(G(x  +  s)  —  G))  (5) 

where  G  is  the  average  value  of  the  field.  Hence,  analyzing  the  field  using  the  SF  only  might 
not  reveal  its  scaling  properties  if  P  <  1.  Yaglom  (1987,  p  390)  shows  that  the  spectrum  is  the 
Fourier  transform  of  R(s)  for  P  <  1  and  the  Fourier  transform  of  the  SF  for  P  >  1.  For  this  reason, 
it  is  expedient  to  verify  the  scaling  of  a  field  using  the  spectrum  first.  Then  depending  on  the 
value  of  P,  one  could  use  either  the  structure  function  (for  /?>  1  )  or  singularity  analysis 
(Davis  et  al.,  1996;  Boufadel  et  al.,  2000)  (for  p  <  1)  to  quantify  multiffactality. 

The  discussion  above  is  general  and  does  not  depend  on  the  adopted  multifractal  model. 
However,  the  adoption  of  such  a  model  simplifies  the  analysis  and  allows  generations  of 
stochastic  realizations.  We  will  use  for  this  purpose  the  universal  multifractal  model  developed 
by  Schertzer  and  Lovejoy  (1987).  The  Schertzer  and  Lovejoy  (S&L)  model  is  a  generalization  of 
the  log  normal  model  (Kolmogorov,  1962;  Obukhov,  1962;  Yaglom,  1966;  Monin  and  Yaglom, 
1975,  p  614).  It  relies  on  using  the  a-stable  family  distribution,  which  is  a  generalization  of  the 
Gaussian  distribution,  to  include,  for  example,  infinite  variance  distributions.  The  a-stable 
distribution  depends  on  four  parameters  (Grigoriu  1995;  Nikias  and  Shao  1995;  Uchaikin  and 
Zolotorev  1999),  but  in  multifractal  studies,  two  of  the  parameters  are  fixed  (Boufadel  et  al. 
2000).  Hence,  only  two  parameters  are  of  interest:  a  and  Ci.  The  parameter  a  is  such 
0  <a<  2  ,  the  upper  limit  being  the  Gaussian  distribution.  As  a  decreases,  the  frequency  of 
sudden  large  jumps  in  the  random  field  increases.  The  parameter  Ci  is  known  as  the  scale 
parameter.  It  is  a  measure  of  the  width  of  the  distribution,  and  is  equal  to  half  the  variance  when 
a  =  2.  As  Ci  increases,  the  magnitude  of  the  sudden  large  jumps  increases  (Boufadel  et  al.  2000). 
Studies  have  shown  that  the  S&L  model  is  flexible  in  fitting  observed  data  (Lavallee,  1991; 
Schmitt  et  al.  1995;  Liu  and  Molz,  1997;  Tennekoon  et  al.,  2003).  In  the  S&L  model  the  structure 
function  exponent  (Eq.  3)  is  given  by: 

§(?)=  qH-  —  -  (qa-q)  (6) 

(a-1) 

where  H  is  a  parameter  that  is  commonly  in  the  range  [0,1].  Hence,  knowledge  of  the  values  of 
the  parameters  a,  Cl5  and  H  is  sufficient  to  characterize  a  field  and  to  generate  stochastic 
realizations  of  the  field. 

One  simple  way  to  generate  multiffactals  is  by  multiplicative  cascades.  This  approach  was 
originally  developed  to  model  turbulent  flow,  where  energy  from  large  eddies  is  progressively 
dissipated  into  smaller  and  smaller  eddies.  Figure  6  illustrated  the  process.  Starting  with  a  field 
that  is  homogeneous  over  the  entire  region  interest,  the  value  of  the  field  (magnetic  susceptibility 
in  our  case)  is  statistically  redistributed  to  finer  scales  according  to  probabilities  governed  by  the 
selected  multiffactal  parameters.  Because  the  total  volume  under  the  curve  is  conserved, 
increasingly  high  peaks  develop,  separated  by  regions  of  small  values  causing  the  field  to  appear 
sparse  at  some  locations.  By  judicious  selection  of  model  parameters,  a  wide  variety  of  textures 
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can  be  simulated.  Figure  6  shows  multifractal  simulation  by 
discrete  cascade,  but  it  is  also  possible  to  synthesize 
continuous  multifractals  using  Fourier  methods  (Wilson  et  al, 
1991;  Pecknold,  1993;  Tchiguirinskaia  et  al.,  2000),  which  is 
the  method  we  employed  in  our  simulations. 

Whether  simulated  by  discrete  or  continuous  cascade,  the 
key  is  to  first  use  the  data  to  obtain  the  multifractal 
parameters  a,  Ci,  and  H.  We  will  illustrate  how  this  is  done 
our  discussion  of  our  experiments  with  synthetic  data  sets. 
But  first  we  discuss  the  justification  for  using  multifractals  to 
model  magnetic  data. 

Multifractal  Modeling  and  Magnetic  Data 

Our  choice  of  fractals  as  a  model  for  soil  and  rock  properties 
is  motivated  by  previous  works  in  the  field.  Pilkington  and 
Todoeschuck  (1993)  analyzed  (in  the  horizontal  direction) 
magnetic  data  in  Canada  using  Fourier  techniques,  and  found 
the  spatial  power  spectrum  to  decay  as  k"p,  where  p  is  a  real 
number  such  that  2  <  p  <  4.  Such  a  power  law  behavior 
evidences  scaling  regime.  These  findings  were  further 
confirmed  by  the  authors  in  another  publication  (Pilkington 
and  Todoeschuck  1995).  Maus  et  al.  (1999)  conducted 
variogram  analysis  on  magnetic  data  in  Canada  and  found 
the  data  to  exhibit  scaling  behavior  as  function  of  the 
horizontal  distance.  Studies  in  the  soil  science  field  have 
shown  that  the  soil  pore  structure  is  fractal  (Perrier  et  al. 
1996;  Baveye  et  al.  1998).  Scaling  behavior  has  also  been 
noted  for  subsurface  properties  (namely  the  intrinsic 
permeability,  K)  as  a  function  of  space.  (Hewett,  1986; 
Neuman,  1990;  Molz  and  Boman,  1993,  1995;  Painter  and 
Paterson,  1994;  Painter  1996,  1998;  Molz  et  al.  1997; 
Boufadel  et  al.,  2000).  Fractality  has  been  found  in  other 
geophysical  fields,  such  as  river  flows  (Gupta  and  Waymire 
1990),  water  in  the  atmosphere  (Davis  et  al.,  1996),  carbon 
content  in  ice  cores  (Schmitt  et  al.  1995),  and  land 
topography  (Lavallee  et  al.,  1993;  Lovejoy  et  al.,  1995). 


Lovejoy  et  al.  (2001)  and  Pecknold  et  al.  (2001) 
published  a  pair  of  papers  investigating  the  multiffactal 
modeling  of  geomagnetic  data  collected  over  the 
Canadian  shield.  Rather  than  model  the  magnetic  field 
directly,  they  chose  to  model  the  underlying  magnetic 
susceptibilities  and  use  their  susceptibility  models  to 


Figure  6.  A  multifractal  cascade  starts 
with  a  uniform  quantity  that  is 
stochastically  subdivided  to  ever  finer 
scales  while  maintaining  the  ensemble 
average.  Eventually  singularities 
(spikes)  and  relatively  quiet  areas 
develop. 
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generate  simulations  of  the  earth's  magnetic  field.  To  do  this  they  made  a  number  of 
assumptions:  (1)  the  magnetization  of  the  geologic  material  is  purely  induced,  no  permanent 
magnetic  moments;  (2)  the  field  can  be  viewed  as  arising  from  a  stratified  earth;  (3)  that  the 
earth's  field  can  be  approximated  as  vertical  (90  degree  inclination,  a  reasonable  approximation 
for  Canada). 

For  this  project  we  also  elected  to  model  the  magnetic  susceptibilities  rather  than  the  magnetic 
field.  Because  we  are  not  interested  in  continental  scale  models,  we  ignored  the  deep  structures 
that  contribute  mainly  to  extremely  large-scale  heterogeneities  in  the  magnetic  field,  and  we 
approximated  the  near-surface  soil  and  bedrock  geology  as  a  single  horizontal  layer  with  a 
multifractal  susceptibility  distribution. 

Experiments  with  Synthetic  Magnetic  Data 

Before  applying  universal  multifractals  to  the  modeling  of  magnetic  data  collected  in  the  field, 
we  tested  whether  we  could  create  a  multifractal  magnetic  field  and  then  analyze  the  synthetic 
data  to  recover  the  multifractal  parameters  uses  to  generate  the  synthetic  data.  The  following 
example  illustrates  the  process  for  a  grid  representing  an  area  half  a  kilometer  on  a  side  with  a 
0.5  m  grid  spacing. 

Generation  of  a  Multifractal  Susceptibilities 

The  generation  of  multifractal  field  is  described  in  detail  by  Tchiguirinskaia  et  al.  (2000), 
Pecknold  et  al.,  (1993)  and  Wilson  et  al.  (1991).  We  summarize  the  steps  and  refer  the  reader  to 
these  publications  for  more  detail. 

Multifractals  are  generated  by  a  multiplicative  cascade,  but  the  same  process  can  be  achieved  by 
adding  quantities  and  then  exponentiating  them  at  the  end.  The  additive  process  is  called  the 
generator,  and  for  this  we  used  random,  negatively- skewed,  Levy  noise  with  the  desired  a  and  Ci 
values.  The  noise  is  then  Fourier  transformed  and  fractionally  integrated  by  mutiplying  by 
|k|_d,“  in  Fourier  space,  where  1/ «  + 1/ «  '=  1  .  After  taking  the  inverse  Fourier  transform 
to  return  to  real  space  the  field  is  exponentiated  to  convert  the  additive  process  to  a  multiplicative 
process.  The  data  are  then  normalized  to  create  a  conservative  multifractal  field.  The  data  are 
then  Fourier  transformed  again  and  fractionally  integrated  by  \k\~K  to  obtain  the  desired 
universal  multifractal  with  parameters  a,  Ci,  and  H.  This  multifractal  field  is  then  scaled  to  give 
the  desired  range  of  values  for  the  magnetic  susceptibilities.  Changing  the  universal  multifractal 
parameters  adjusts  the  smoothness,  and  intermittency  of  the  resulting  multifractal,  making  it 
possible  to  capture  a  wide  range  of  behaviors  while  maintaining  the  fundamental  scaling 
behavior  (the  statistical  nature  of  the  field  is  the  same  at  every  scale).  The  process  is  illustrated 
for  a  1-D  example  in  Figure  7. 

Simulation  the  Magnetic  Field 

To  go  from  the  multifractal  model  of  magnetization  to  the  magnetic  field,  we  used  a  simple 
horizontal  sheet  model.  This  is  reasonable  given  that  it  is  the  geologic  variability  of  the  shallow 
soil  and  bedrock  that  produces  the  short-wavelength  anomalies  can  obscure  the  signature  of 
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Figure  7:  The  figure  illustrate  the  creation  of  a  multifractal  using  Fourier  methods  for  the  1-D 
case.  From  top  to  bottom:  Levy  noise  is  generated  with  the  desired  a  and  Ci  values.  The  noise  is 
converted  to  an  additive  process  by  partial  integration  in  the  Fourier  domain.  The  additive 
process  is  converted  to  a  multiplicative  process  by  exponentiation  after  returning  to  the  spatial 
domain.  Finally,  a  second  partial  integration  in  the  Fourier  domain  yields  a  multifractal  with 
the  desired  value  ofH.  (a  =  1.8,  Ci  =  0.5,  H  =  0.4) 


Figure  8:  Model  for  a  0.5  m-thick  horizontal  layer  with  1-m  horizontal  cells  having  a 
multifractal  magnetic  susceptibility  distribution.  The  susceptibility  values  here  are  scaled 
between  0-1,  but  can  be  adjusted  to  have  any  desired  range. 


13 


small  UXO  targets.  For  a  sheet  model  the  magnetization  varies  only  in  the  horizontal  direction 
(Figure  8).  Assuming  purely  induced  magnetization,  the  magnetic  field  at  any  elevation  above 
ground  level  can  be  calculated  in  the  Fourier  domain  using  the  equation  (Blakely,  1996): 

F[AT}=  F\M}{2nCm0m0fe(lkUo){e 


idkko)  (  „~\k\Zl  _ e~\k\z2 


(7) 


z0<zl,zl<z 


2  ’ 


where  F[AT\  is  the  Fourier  transform  of  the  magnetic  anomaly,  F[M]  is  the  Fourier 
transform  of  the  magnetization,  Cm  is  a  constant  (10'7  henry/m),  0m  and  0f  are  the 
orientation  vectors  for  the  magnetization  and  the  earth's  field,  respectively,  k  is  the  wavenumber, 
z0,  is  the  height  above  the  surface  of  the  measurements,  and  zu  z2,  are  the  depths  below  the 
surface  to  the  top  and  bottom,  respectively,  of  the  magnetic  layer.  The  multifractal 
magnetizations  can  then  be  scaled  to  match  the  range  of  amplitudes  encountered  in  the  field. 
Given  a  multifractal  realization  of  the  magnetic  susceptibility,  the  corresponding  magnetic  field 
can  be  calculated  at  any  altitude.  Figure  9a  shows  a  simulation  of  magnetic  susceptibilities  over  a 
square  area  0.5  km  on  a  side  with  a  grid  spacing  of  0.5  m.  Figure  9b  shows  the  resulting 
magnetic  field  calculated  at  a  height  of  2.0  m  assuming  an  a  multifractal  layer  of  magnetic 
susceptibilities  on  a  grid  spacing  of  0.5  m  created  using  universal  multifractal  parameters  H=0.3, 
a  =  1.8,  and  Ci  =  0.050,  and  assuming  an  inducing  magnetic  field  with  an  inclination  of  62°,  and 
a  declination  of  10°.  To  apply  our  methodology  to  field  data,  we  must  first  demonstrate  that  we 
can  successfully  recover  the  original  multifractal  model  parameters  from  the  magnetic  data  for 
this  synthetic  example. 


Figure  9:  (A)  Multifractal  distribution  of  magnetic  susceptibilities  (S.I.)  simulated  on  a  0.5  m 
grid.  (B)  The  corresponding  total  magnetic  field  calculated  at  2  m  above  the  surface  (nT). 


Recovery  of  the  Universal  Multifractal  Parameters 

We  have  shown  how  to  produce  multifractal  distributions  of  magnetic  susceptibilities  and  then 
use  these  to  generate  a  magnetic  field,  but  given  a  magnetic  field,  how  do  we  recover  the 
parameterization  of  of  the  underlying  multifractal  distribution?  This  is  the  problem  we  face  with 
field  data,  so  we  need  to  verify  that  our  methodology  performs  correctly  for  synthetic  data. 

Starting  with  the  simulated  magnetic  data  (Figure  9b),  the  first  step  is  to  test  that  a  multifractal 
model  is  appropriate  by  checking  that  the  data  are  scaling.  This  we  do  by  checking  that  a  log-log 
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plot  of  the  Fourier  spectrum  of  the  data  are  linear.  A 
plot  of  the  spectrum  of  the  magnetic  susceptibilities  is 
linear  (Figure  10),  but  the  Fourier  spectrum  of  the 
calculated  magnetic  field  clearly  is  not  (Figure  11). 

This  is  because  the  field  is  calculated  at  a  height  of  2  m 
above  ground  level  to  simulate  data  collection  by 
helicopter,  and  calculation  of  the  magnetic  field  at 
altitudes  above  the  source  is  equivalent  to  upward 
continuation,  which  is  accomplished  in  the  Fourier 
domain  by  multiplying  by  a  factor  of  e~^ '  , 
rendering  the  spectrum  non-linear  by  preferentially 
attenuating  the  higher  frequencies.  The  obvious 
remedy  is  to  downward  continue  the  data  back  to 
ground  level  by  multiplying  by  e^z  in  the  Fourier 
domain.  Upward  continuation  is  stable,  because 
multiplication  by  a  negative  exponential  is  a 
smoothing  operation.  Downward  continuation, 
however,  involves  multiplication  by  a  positive 
exponential,  and  is  notoriously  unstable.  Any  noise 
present  in  the  data,  even  numerical  noise,  is 
magnified  exponentially  by  the  downward 
continuation  filter.  Traditionally,  the  way  around  this 
is  to  first  low-pass  filter  the  magnetic  data  before 
downward  continuing  it,  but  low-pass  filtering  will 
strongly  alter  the  underlying  spectrum. 

A  recent  paper  by  Xu  et  al.  (2007)  proposed  a  new 
approach  to  downward  continuation  that  is 
numerically  stable,  and  dispenses  with  low-pass 
filtering  at  the  expense  of  additional  computational 
time.  The  idea  is  straightforward;  the  magnetic  field 
at  the  measurement  elevation  is  used  as  an  initial 
guess  for  the  magnetic  field  at  the  ground  level,  then 
upward  continued  (stable  operation)  and  compared 
with  the  measured  data.  The  field  at  the  ground  level 
is  then  iteratively  adjusted  until  the  upward  continued 
field  and  the  measured  data  agree  within  an 
acceptable  tolerance.  Using  this  approach  to 
downward  continuation,  no  spurious  high  frequencies 
are  introduced.  Figure  12  shows  the  spectrum  of  the 
downward  continued  magnetic  data,  which  is  clearly 
linear  with  the  same  slope  as  the  spectrum  of  the 
simulated  magnetic  susceptibilities,  except  at  the 
highest  frequencies  where  information  was  lost 
through  the  original  upward  continuation  operation. 


V 
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Figure  10:  Fourier  power  spectrum  of 
the  simulated  magnetic  susceptibilities. 


Averaged  Power  Spectral  Density  (X-Direction) 


Figure  11:  Fourier  power  spectrum  of 
magnetic  field  at  a  flight  height  of  2  m 

Averaged  Power  Spectral  Density  (X-Direction) 


Figure  12:  Fourier  spectrum  of  magnetic 
data  downward  continued  to  the  surface, 
above  the  surface. 
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The  Fourier  spectrum  of  the  downward  continued  magnetic  data  shows  that  the  underlying 
magnetic  susceptibilities  are  scaling,  therefore  a  fractal  model  is  appropriate.  Thus,  the  first  step 
in  analyzing  field  data  will  be  to  downward  continue  the  data  and  check  the  Fourier  spectrum. 

After  determining  the  data  are  scaling,  the  next  step  is  to  find  the  corresponding  universal 
multifractal  parameters  a,  Ci,  and  H.  For  a  conservative  field  (H=0),  the  structure  function  (Eq. 
6)  becomes: 


K(q)=  (qa~q).  (8) 

The  parameters  a,  Ci  can  be  estimated  using  the  double  trace  moment  method  (DTM)  (Lavallee 
et  al,  1992,  1993;  Boufadel  et  al.,  2000).  In  this  technique  the  data  are  first  raised  to  a  power  r|, 
then  averaged  over  a  scale  X  (the  term  “dressed”  is  often  used  in  the  literature).  The  resulting 
field  will  have  the  moment  scaling  function,  K(q,r|)  for  a  universal  multifractal  given  by: 

K{q,q)=  qaK{q)=  qa-^1—(qa-q),  a*  1.  (9) 

(a-1) 


We  obtain  an  estimate  of  K(q,r|)  from  the  data  using  Eq.  3  to  calculate  the  qth  moment,  for  the 
data  averaged  using  different  values  of  r|.  The  slope  of  a  log-log  plot  of  K(q,rp  versus  r\  yields 
an  estimate  of  the  parameter  a.  Using  different  values  of  q  should  yield  parallel  lines  with  the 
same  slope,  serving  as  a  check  on  the  a  estimate.  Given  a,  we  can  deduce  Ci  from  same  plot 
from  the  value  of  K(q,q)  at  r\  =  1  and  Eq.  9.  There  is  one  caveat.  The  process  breaks  down  in 
practice  for  higher  order  moments  (typically,  q  >  4)  where  one  or  two  high  data  values  dominate 
the  ensemble  averages,  and  can  also  break  down  for  extremely  low  values  of  q  due  to  noise  or 
discretization  error  (Pecknold  et  al.,  2001).  Avoiding  these  extremes,  the  DTM  technique  yields 
good  estimates  of  a  and  Ci. 

For  a  non-conservative  field  (H  >  0),  one  has  a  third  parameter  to  estimate,  H.  The  DTM 
technique  applied  to  non-conservative  field  yields  incorrect  values  for  a  and  Ci  If  the  value  of  H 
were  known,  one  could  convert  the  non-conservative  field  into  a  conservative  one  by  fractional 
differentiation  of  order  H  (equivalent  to  multiplying  by  kH  in  Fourier  space,  where  k  is  the 
wavenumber).  Not  knowing  the  value  of  H,  we  employ  a  “bootstrap  method.”  We  take  the 
increment  of  the  data  (which  is  equivalent  to  fractionally  differentiating  by  assuming  H  =  1),  and 
use  the  DTM  to  estimate  a  and  Ci.  The  DTM  technique  will  yield  reasonable  estimates  as  long 
as  H  <  1 .  For  a  conservative  field,  the  slope,  pcon ,  of  a  line  fitted  to  a  log-log  plot  of  the  Fourier 
power  spectrum  is: 


/U=  l-*(2). 


(10) 


where  K(q)  is  given  by  Eq.  9.  Using  the  values  for  a  and  Ci  obtained  from  the  DTM  results  for 
the  differenced  data,  we  then  calculate  pconand  compare  with  the  P  value  obtained  from  the 
power  spectrum  of  the  original  data.  The  value  of  H  can  then  be  estimated  as: 
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(11) 


This  value  of  H  is  then  used  to  fractionally  differentiate  the  original  data,  the  DTM  technique 
applied  again  to  obtain  new  estimates  of  a  and  Ci ,  a  new  value  of  H  obtained  from  Eq.  11,  and 
the  whole  process  repeated  until  the  parameter  estimates  converge  (Pecknold  et  al.,  2001), 
usually  within  a  few  iterations.  Applying  this  method  to  the  downward  continued  magnetic  data 
in  our  example,  we  obtain  parameter  values  of  H=0.4,  a  =  1.7,  and  Ci  =  0.059,  which  are  quite 
close  to  the  parameters  used  to  generate  the  multifractal  magnetic  susceptibilities  (H=0.3,  a  = 

1.8,  and  Ci  =  0.050).  Thus,  we  have  described  a  procedure  for  fitting  a  universal  multifractal 
model  to  total  field  airborne  magnetic  data  and  generating  new  realizations  with  the  same 
statistical  behavior  at  scales  both  smaller  and  larger  than  the  original  data  set.  We  will  now  apply 
the  methodology  to  several  field  data  examples. 

Case  Histories  of  Multifractal  Simulation  of  Field  Data 

We  evaluated  our  methodology  using  three  different  test  sites  representing  low,  medium,  and 
high  levels  of  geologic  noise  relative  to  the  magnetic  signal  strength  of  buried  UXO.  For  the 
low-noise  case  we  selected  the  Sierra  Army  Depot  in  California.  For  the  medium-noise  site  we 
selected  Fort  Ord,  California,  and  for  the  high-noise  site  we  selected  Isleta  Pueblo,  New  Mexico. 
We  will  present  the  case  histories  in  order  of  a  increasing  geologic  noise. 

Case  History:  Sierra  Army  Depot,  CA  -  Minimal  Geologic  Noise 


The  Sierra  Army  Depot  is  located  in  southern 
California,  just  across  the  border  from  Reno, 
Nevada.  A  team  from  Oak  Ridge  National 
Laboratory  conducted  an  airborne  geophysical 
survey  in  January  and  February  of  2003, 
surveying  the  Honey  Lake  area  (3466  acres),  and 
the  East  Shore  (1168  acres)  using  the  helicopter- 
based  ORAGS™  Arrowhead  system  (Figure  13) 
which  is  comprised  of  a  linear  array  of  eight 
Scintrex™  Cesium  vapor  magnetometers  (ORNL, 
2003)  spacing  1.75  m  apart  such  that  each  pass 
cuts  a  swath  just  over  12  m  wide.  The  surveys 
were  flown  at  an  average  flight  height  of  1.75  m 


above  ground  level  for  Honey  Lake,  and  2.21m  Figure  13:  The  ORNL  Arrowhead  system  in 
for  the  East  Shore.  flight  at  SIERRA  Army  Depot. 


According  to  California  Groundwater  Bulletin  118  (California,  2003),  the  bedrock  under  the 
Honey  Valley  Basin  consists  of  Plio-Pleistocene  and  Pleistocene  age  basalt  flow,  which  are  likely 
to  be  highly  magnetic,  but  these  flows  are  overlain  by  lake  and  near-shore  sediment  deposits  up 
to  213  m  (700  ft)  thick.  Comprised  of  sand,  clay,  and  silt,  these  non- magnetic  overlying 
sediments  provide  a  large  separation  between  near-surface  UXO  targets  and  magnetic  bedrock. 
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Sierra  Army  Depot:  Honey  Lake  and  East  Shore 

Residual  Magnetics 

ORNL  Helicopter  Survey,  January  -  February  2003 
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GEOPHYSICS 


Figure  14:  Residual  magnetic  data  from  SIERRA  Army  Depot.  Notice  the  scale.  The 
background  magnetic  variation  away  from  the  UXO  is  only  a  few  nT  (from  ORNL,  2003). 
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The  result  is  a  site  with  virtually  no  magnetic  noise  originating  from  near- surface  geologic 
sources.  Figure  14  shows  the  residual  total  magnetic  field  after  removing  the  regional  N-S  trend. 
Notice  that  the  background  away  from  the  UXO  anomalies  varies  by  at  most  a  few  nanotesla. 
This  level  of  compensation  for  helicopter-related  noise  sources  (rotor  noise,  flight  level 
fluctuations,  GPS  time  lag,  tilt-roll-yaw,  heading  error,  etc.)  is  far  easier  to  achieve  for  residual 
magnetic  data  then  for  total  field.  High-pass  filtering  down  the  lines  followed  by  adjustments  to 
bring  each  trace  to  zero  mean  eliminates  most  of  the  compensation  error.  Unfortunately,  it  also 
distorts  the  frequency  content  of  any  remaining  geologic  noise,  making  it  impossible  to  apply  the 
multifractal  methodology. 

We  analyzed  a  subset  of  the  ORNL  magnetic  data,  selecting  a  region  with  as  few  UXO 
anomalies  as  possible  to  characterize  the  geologic  background  without  interference,  and 
attempted  to  reprocess  the  total  field  data  to  recover  the  geologic  background.  Figure  15  shows 
the  total  field  data  for  our  1-km  test  area  is  dominated  by  a  N-S  trend  created  by  deep-seated 
magnetic  sources.  After  removing  the  regional  N-S  trend  by  fitting  and  subtracting  a  second 
order  polynomial  surface  the  magnetic  “background”  was  dominated  by  small  amounts  of  line- 
to-line  compensation  error  (rotor,  elevation  fluctuations,  etc.),  which  masks  any  anomalies 
created  shallow  magnetic  rock  or  soil.  It  is  possible  that  geologic  noise  of  less  roughly  10-15nT 
is  present,  but  further  reducing  the  compensation  error  would  require  microleveling  the  flight 
lines.  This  would  require  data  from  E-W  tie  lines,  which  were  not  collected. 

We  conclude  that  at  sites  such  as  SIERRA  Army  Depot,  where  the  geologic  background  is 
minimal,  further  improvement  in  magnetic  detection  and  discrimination  will  come  with  reduction 
in  system  noise,  not  synthetic  modeling  of  the  magnetic  signature  of  the  geology.  As  we  shall 
see,  this  is  not  the  case  for  our  two  remaining  case  histories. 


Figure  15:  (Left)  The  total  magnetic  field  for  a  1-km  test  portion  of  the  Sierra  Army  Depot 
magnetic  data  is  dominated  by  a  N-S  trend  created  by  deep-seated  geologic  sources.  (Right) 
With  the  trend  removed  line-to-line  residual  compensation  error  dominates.  Shallow  geologic 
sources  are  too  weak  to  identify.  Color  scale  is  in  nT. 


19 


Case  History:  Fort  Ord,  CA  -  Moderate  Geologic  Noise 


Fort  Ord  is  situated  near  Monterey  Bay, 

California,  roughly  130  km  (80  miles)  from 
San  Francisco,  and  close  to  the  Pacific  coast. 

The  base  was  officially  closed  in  September, 

1994,  however  remediation  activities  continue 
at  the  site,  including  the  identification  and 
removal  of  UXO.  In  2005,  the  ORNL  team 
conducted  both  helicopter  magnetic  and 
electromagnetic  surveys  of  approximately 
1,281  hectares  (3166  acres)  of  Fort  Ord 
looking  for  clusters  of  anomalies  associated 
with  targets  and  ranges  (ORNL,  2005). 

The  geology  is  comprised  of  mesozoic  granite 
and  metamorphic  rocks  overlain  by  younger 
alluvial  fan,  lake,  and  fluvial  deposits  (Fort 
Ord  Primer,  2008).  Consequently  we  would 
expect  the  bedrock  to  be  magnetic,  but  the 
overlying  sediments  to  be  non-magnetic,  as  we 

saw  was  the  case  for  the  Sierra  Army  Depot.  Magnetic  susceptibility  measurements  the  ORNL 
team  on  grab  samples  of  the  surface  sediments  confirmed  that  the  overlying  sands  are  non¬ 
magnetic.  And,  as  at  are  Sierra,  the  Fort  Ord  magnetic  data  are  dominated  by  a  N-S  trend 
probably  associated  with  the  underlying  bedrock.  Apart  from  some  magnetic  background 
associated  with  debris,  the  researchers  noted,  “Numerous  other  moderate  amplitude  responses 
exist  within  the  survey  area,  but  these  are  likely  more  geologic  in  origin.”  (ORNL,  2005)  These 
sinuous  anomalies  appear  to  coincide  with  Fort  Ord's  ridges  and  trough  topography  (Figure  16). 
ORNL  researchers  speculated  that  magnetic  sediments  may  have  washed  down  from  the  nearby 
mountain  ranges  can  accumulated  in  the  troughs  (personal  communication,  2008),  although 
confirmation  of  this  hypothesis  would  require  ground  follow-up. 

Multifractal  analysis  of  the  Fort  Ord  magnetic  data  was  complicated  by  two  factors  in  addition  to 
the  aforementioned  anisotropy.  First,  vegetation  coupled  with  the  uneven  terrain  made  flying 
close  to  the  ground  more  difficult  here  than  at  Sierra  Army  Depot.  The  average  flight  height  at 
Sierra  was  less  than  2.0  m  above  ground  level;  at  Fort  Ord  it  was  3.5  m.  The  magnetic  response 
from  a  small,  localized  source  falls  off  with  the  cube  of  distance,  so  nearly  doubling  the  flight 
height  means  a  factor  of  eight  loss  in  sensitivity.  More  importantly  for  the  multifractal  analysis, 
it  means  an  increase  in  the  smallest  scale  that  can  be  resolved.  The  second  challenge  in 
analyzing  these  data  was  the  survey  design.  To  economize,  the  ORNL  team  flew  every  other 
line.  The  Arrowhead  system  (Figure  16)  records  data  at  a  rate  of  120  data  points  per  second  from 
eight  magnetometers  spaced  1.7  m  apart  perpendicular  to  the  flight  direction.  Flying  over  Ford 
Ord  at  an  average  speed  of  roughly  20  m/s,  resulted  in  an  extremely  fine  down-the-line  sampling 
interval  of  less  than  20  cm  for  lines  spaced  1.7  m  apart.  But  by  skipping  every  other  line  the 
final  result  was  a  data  set  alternating  eight  tightly  spaced  lines  with  12  m  gaps  between  the  line 
groups  (Figure  17). 


Figure  16:  ORNL  Arrowhead  system  in  operation 
at  Fort  Ord.  Notice  the  ridge  and  trough 
topography,  which  is  apparently  correlated  with 
the  anisotropic  magnetic  background. 
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Figure  17:  Residual  total  magnetic  field  from  Fort  Ord.  Notice  the  strong  N-S  trend  from  deep- 
seated  sources,  the  12-m  gaps  created  by  flying  every  other  line,  and  the  strong,  sinuous 
anisotropy  created  by  near-surface  geologic  sources.  ( From  ORNL,  2005).  Note:  Units  for  this 
map  are  in  State  Plane  feet;  each  grid  cell  is  1000ft  on  a  side.  Consequently,  the  test  area  is 
roughly  3x3  cells  and  is  shown  as  a  shaded  region  above. 

Multifractal  analysis  of  the  Fort  Ord  data  began  with  the  selection  a  1  km2  subset  of  the  data  with 
a  minimum  of  UXO  and  other  anthropogenic  magnetic  anomalies.  Figure  18  shows  the  test  area 
in  more  detail,  scaled  to  emphasize  the  magnetic  background.  Even  in  this  relatively  clean 
potion  of  the  Fort  Ord  data  set  the  few  UXO  items  present  produced  by  far  the  largest  anomalies. 
Fortunately  for  our  goal  of  background  estimation,  the  larger  UXO  items  contributed  to  only  a 
small  portion  of  the  test  area  data,  and  after  clipping  the  UXO  anomaly  peaks  they  contributed 
negligibly  to  the  power  spectrum.  The  contribution  of  smaller  UXO  items  that  buried  in  the 
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Figure  18:  Total  field  magnetic  data  for  1  sq  km  test  area  selected  from  the  Fort  Ord  data  set 
(see  Fig.  17)  as  a  region  with  minimal  UXO.  Although  some  UXO  are  present,  the  anisotropic 
magnetic  background  is  dominant  for  most  of  this  region.  ( Color  scale  is  in  nT) 


background  noise,  however,  is  uncertain.  In  the  analysis  that  follows  we  implicitly  assume  that 
the  sinuous,  anisotropic  anomalies  presumably  associated  with  magnetic  sediments  dominate. 
Note  that  for  this  test  area  the  data  have  been  interpolated  to  a  2-m  grid,  including  interpolation 
across  the  data  gaps  created  by  the  alternating  flight  lines,  which  may  affect  the  information  and 
the  shortest  wavelengths,  but  was  unavoidable  in  this  case. 

The  second  step  in  the  multiffactal  analysis,  after  selecting  a  relatively  clean  subset  of  the  data,  is 
to  downward  continue  the  data  to  the  ground  surface,  then  examine  the  Fourier  power  spectrum 
(Figure  19).  The  result  clearly  shows  fractal  scaling  (linear  decrease  with  frequency)  except  for 
the  highest  frequencies,  which  cannot  be  fully  restored  by  downward  continuation.  Following 
the  methodology  described  earlier,  we  now  calculate  the  double  trace  moment  to  determine  the 
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Averaged  Power  Spectral  Density  (X-Direction) 


Averaged  Power  Spectral  Density  (Y-Direction) 


Figure  19:  Power  spectrum  for  the  1  sq.  km 
subset  of  the  Fort  Ord  magnetic  data.  The 
linear  decrease  of  the  logarithm  of  the  spectrum 
with  the  logarithm  of  frequency  is  evidence  of 
fractal  scaling.  Linearity  breaks  down  at  the 
highest  frequencies,  which  cannot  be  accurately 
restored  by  downward  continuation. 


Figure  20:  The  double  trace  moment  method 
applied  to  the  data  for  various  moments,  q, 
recovers  two  of  three  required  universal 
multifractal  parameters:  alpha  and  cl. 


universal  multifractal  parameters  (Figure  20).  The  fits  show  the  expected  departure  from 
linearity  (divergence  of  the  moments)  for  high  values  of  r|  (from  Eq.  9).  The  parallel  linear 
slopes  for  different  values  of  q,  however,  lends  confidence  to  the  obtained  values  of  a  =  1.8  and 
cl  =  0.08.  The  value  of  H  obtained  as  described  in  the  methods  section  was  H=0.3. 

We  now  have  the  parameters  required  to  simulate  the  magnetic  background.  However,  the 
simulation  methodology  assumes  isotropy,  which  is  clearly  a  poor  assumption  for  this  site.  We 
have  recently  begun  exploring  methods  for  incorporating  anisotropy  in  our  simulations,  though 
this  goes  beyond  the  original  scope  of  this  pilot  project.  One  straightforward  though  not  very 
general  approach  to  simulating  anisotropy  is  to  use  different  values  for  Hx  and  Hy  in  the  final 
fractional  integration  that  converts  the  conservative  multifractal  into  an  non-conservative 
multifractal  (means  varies  with  scale).  Figure  21a  shows  the  magnetic  data  from  our  test  area 
compared  with  one  realization  of  an  isotropic  simulation  (Figure  21b)  using  the  fitted  universal 
multifractal  parameters,  and  two  anisotropic  (Figure  21c,d)  realizations  using  anisotropic 
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Figure  21:  (a)  Subset  of  the  Fort  Ord  total  field  magnetic  data;  (b)  Isotropic  multifractal 
simulation  using  the  fitted  parameters  (see  text);  (c  &  d)  Anisotropic  multifractal  simulations. 
Color  scale  is  in  nT. 


fractional  integration,  with  the  parameters  selected  by  eye  to  achieve  a  reasonable  match. 

We  conclude  from  this  case  history  that  the  magnetic  background  at  Fort  Ord  is  anisotropic 
multifractal,  and  reasonable  simulations  can  be  achieved  with  a  simple  modification  to  the 
simulation  of  isotropic  universal  multifractals.  However,  the  ad  hoc  nature  of  this  anisotropic 
simulation  is  unsatisfying.  A  more  general  approach  to  the  characterization  and  simulation  of 
anisotropic  universal  multifractals  is  possible  (e.g.,  Cheng,  2004;  Lewis,  1999;  Lovejoy  et  al., 
2001),  which  can  incorporate  both  anisotropic  stratification  as  well  as  change  in  the  direction  of 
anisotropy  (rotation)  with  scale,  and  this  will  be  the  subject  of  our  future  investigations. 
Nonetheless,  our  simple  anisotropic  model  has  produced  quite  reasonable  simulations  of  the 
magnetic  background  at  Fort  Ord. 
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Case  History:  Pueblo  of  Isleta,  NM  -  Strong  Geologic  Noise 

The  Pueblo  of  Isleta  is  located  near  Albuquerque, 

New  Mexico.  Established  in  the  1300's,  the  Native 
American  settlement  comprises  now  Oraibi  and 
Chicale.  Over  the  years  portions  of  this  land  have 
been  contaminated  with  UXO  by  DoD  training 
missions.  Researchers  from  both  ORNL  and  the 
Naval  Research  Laboratory  have  conducted 
separate  demonstrations  wide-area  assessment  at 
the  S 1  target  area  using  of  ultra-low  airborne 
magnetic  mapping.  From  the  point  of  view  of 
flight  logistics,  the  site  is  ideal  -  flat  with  minimal 
vegetation  (Figure  22)  -  but  still  represents  a 
difficult  challenge  for  magnetic  detection  and 
discrimination  of  UXO  because  of  the  strong  Figure  22:  The  MTADS  airborne  system 
magnetic  background.  The  surface  material  is  thin  deployed  at  Pueblo  of  Isleta,  NM.  The  flat 
veneer  of  Holocene  and  Pleistocene  alluvial  and  terrain  and  minimal  vegetation  made  it 
colluvial  deposits,  which  are  probably  only  possible  to  survey  this  site  with  average  flight 

moderately  magnetic,  but  are  underlain  by  a  series  heights  below  2  m  ( Nelson  et  al.,  2003). 
of  ancient  basalt  flows  which  are  highly  magnetic. 

Quoting  from  an  ORNL  report  (ORNL,  2003): 

“39%  of  the  ORAGS  detections  that  were  dug  in  S-01  were  classed  as  ‘no  finds.’  This 
value  is  high  in  comparison  to  other  surveys  where  more  standard  field  excavation 
techniques  were  used  (for  example,  <3%  at  BBR,  Van  et  al.,  2004).  This  number  is 
artificially  high  in  part  because  excavation  radii  did  not  go  beyond  lm.  One  meter  is 
approximately  the  average  separation  between  the  excavated  items  and  their  predicted 
location  based  on  ORAGS-Arrowhead  data.  Other  factors,  including  localized  zones  of 
rock  or  soil  with  high  magnetic  susceptibility  (‘hot  rock/dirt’)  may  have  contributed  to 
the  high  rate  of  ‘no  finds’  in  area  S-01.” 

Similarly,  the  NRL  survey  report  (Nelson  et  al.,  2003)  reports  in  their  results  table  instances  of 
“target  lost  in  the  geology.”  Thus,  magnetic  noise  was  a  major  problem  for  these  surveys. 

The  sheer  density  of  UXO  items  at  the  SI  site  also  presented  a  challenge  for  our  multifractal 
analysis,  making  it  difficult  to  distinguish  the  geologic  background  (Figure  23).  Again,  we 
selected  a  portion  of  the  dataset  away  from  the  main  target  area,  but  even  for  this  subset,  the 
magnetic  signatures  of  the  larger  UXO  present  dominated  the  data  (Figure  24a).  Consequently, 
we  used  UXOLab  (UBC,  2008)  to  model  and  remove  the  effect  of  the  35  largest  dipole 
anomalies,  clearly  too  large  and  too  localized  to  be  of  geologic  origin.  Our  multifractal 
characterization  of  the  residual  assumes  that  the  remaining  magnetic  field  is  dominated  by  the 
geologic  noise,  even  though  many  smaller  UXO  items  remain  (Figure  24b).  Ideally  for  our 
purposes,  the  contractor  would  collect  background  data  in  the  same  geologic  setting  away  from 
UXO-contaminated  sites.  Absent  uncontaminated  background  data,  we  proceeded  with  our 
analysis,  recognizing  that  the  resulting  multifractal  simulation  will  tend  to  overestimate  the 
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Figure  23:  Residual  magnetic  data  for  Pueblo  oflsleta  (ORNL,  2003).  Shaded  box  shows  the 
area  selected  for  multifractal  analysis. 
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Figure  24:  (a)  Surface  plot  of  the  subset  of  the  Isleta  Figure  25:  Power  spectra  of  the  data 
magnetic  data  (boxed  area  in  Fig.  22)  demonstrates  subset  show  scaling  except  beyond  the 
that  even  away  from  the  S-01  target  UXO  dominates  fitting  range  where  downward 
the  total  magnetic  field,  (b)  The  same  data  with  the  continuation  has  affected  the  highest 
35  largest  dipole  anomalies  modeled  and  removed,  frequencies. 

Smaller  UXO  items  remain ,  but  now  the  magnetic  of 

number  of  background  spikes  (largely  affecting  the  cl  parameter).  We  will  return  to  this  point 
this  point  in  our  summary  discussion  at  the  end.  The  Fourier  power  spectra  (Figure  25)  of  this 
subset  of  the  Isleta  data  show  evidence  of  fractal  scaling,  though  the  spectra  depart  from  linear  at 
the  highest  frequencies,  where  downward  continuation  of  the  aerial  magnetic  data  could  properly 
recover  spectral  energies.  We  placed  more  credence  in  the  spectrum  for  the  Y-direction  as  the  X- 
direction  spectrum  also  has  some  contamination  by  line-to-line  leveling  error. 
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0.35 


Figure  26:  Double  trace  moment  method 
applied  to  the  Isleta  test  area  data. 


Figure  27:  Direct  fitting  of  the  structure  function 
(Eq.  6).  The  concave  downward  curvature  is 
evidence  that  the  data  are  multifractal. 


We  then  used  the  DTM  Method  to  obtain  the  values  of  a  =  1.6,  and  cl  =  0.06  (Figure  26).  Then 
fromEq.  10  we  obtained  a  value  of  H  =0.15.  As  a  check,  we  also  fit  the  structure  function  (Eq. 

6)  directly  (Figure  27),  obtaining  the  values  a  =  1.7,  and  cl  =  0.04  and  H  =  0.1.  Differences 
between  the  parameter  values  obtained  using  the  two  approaches  is  too  small  to  affect  the 
multifractal  simulations.  Figure  27  directly  reveals  the  multifractal  nature  of  the  Isleta  data.  For 
a  monofiractal  (ordinary  fractal),  Eq.  6  become  ^(q)  =  qH;  the  structure  function  is  linear.  The 
concave  downward  curvature  in  Figure  27  clearly  shows  that  the  fractal  parameters  are  changing 
with  scale,  hence  the  data  are  multifractal. 

Finally,  Figure  28  shows  the  data  for  our  test  portion  of  the  Isleta  data  and  compared  with  three 
realizations  of  the  multifractal  simulation  using  the  fitted  a,  cl,  and  H  parameter  values.  Recall 
that  the  goal  is  not  to  reproduce  the  exact  noise  pattern  observed  in  the  data,  which  is  only  one 
realization  of  a  complex  geologic  process.  The  goal  is  to  simulate  geologic  noise  with  the  same 
statistical  behavior  as  the  data  across  all  scales.  Figure  28  shows  that  universal  multifractal 
model  captures  this  statistical  behavior,  minus  the  NS  leveling  errors  and  the  remaining  small 
UXO  items.  Furthermore,  once  the  appropriate  multifractal  parameters  have  been  obtained,  we 
are  in  a  position  to  simulate  “Isleta-like”  geologic  noise  on  any  desired  scale,  from  centimeters  to 
hundreds  of  kilometers. 

Discussion  and  Conclusions 

Summary  of  the  Methodology  and  Case  Histories 

The  basic  premise  of  this  work  that  we  can  improve  UXO  detection  and  discrimination  by 
developing  scale-independent  models  for  the  background  noise  created  by  soil  and  bedrock  with 
large  magnetic  susceptibilities.  Quality  Assurance  (QA)  measures  based  on  success  rates  at 
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Figure  28:  Isleta  magnetic  data  ( top  left)  and  three  multifractal  realizations  using  the  universal 
multifractal  parameters  fitted  to  the  test  data. 


UXO  proving  grounds  are  unlikely  to  carry  over  to  field  studies  because  the  local  geology  is 
different,  and  because  geologic  noise  thresholds  change  with  scale.  This  scaling  problem  is 
especially  acute  for  wide-area  surveys  encompassing  10's  or  100's  of  square  kilometers,  hence 
the  need  to  develop  mathematical  models  that  can  generalize  measurements  to  larger  or  smaller 
scales. 

Multifractal  modeling  is  a  well-established  methodology  for  the  representation  of  processes  that 
scale  across  many  orders  of  magnitude.  In  this  pilot  project  we  have  demonstrated  how  to  test 
whether  the  background  magnetic  susceptibilities  are  multiffactal,  how  to  then  obtain  the 
multifractal  parameters  a,  cl,  and  H  to  fit  a  universal  multifractal  model,  and  how  to  use  this 
model  to  generate  simulations  of  the  geologic  background  at  any  scale  desired. 
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We  applied  the  methodology  to  airborne  magnetic  data  sets  collected  at  three  UXO  sites  with 
different  levels  of  geologic  magnetic  noise:  Sierra  Army  Depot,  CA;  Fort  Ord,  CA;  and  Pueblo 
of  Isleta,  NM.  The  geologic  noise  at  Sierra  Army  Depot  was  extremely  low,  less  than  a 
nanoTesla,  well  below  the  magnetic  noise  level  produced  by  the  helicopter  even  after 
compensation.  Improved  UXO  detection  at  such  low-noise  sites  would  require  improved  data 
acquisition  before  multifractal  modeling  would  be  useful. 

Fort  Ord  magnetic  data  showed  higher  levels  of  magnetic  noise,  up  to  +  20  nT  or  more.  Our 
analysis  showed  that  the  data  could  be  fit  with  a  multifractal  model,  but  that  the  magnetic 
background  at  this  site  is  anisotropic.  We  produced  reasonable  anisotropic  simulations  with  a 
fairly  simple  modification  to  our  methodology.  A  more  generalized  approach  would  be  required 
to  incorporate,  for  example,  a  rotational  anisotropy  where  the  direction  of  anisotropy  changes 
with  scale.  Generalized  anisotropy  is  an  issue  we  hope  to  address  in  future  work. 

The  magnetic  data  collected  at  Pueblo  of  Isleta  had  the  highest  levels  of  geologic  background 
noise  of  the  three  data  sets  we  investigated,  with  levels  reaching  more  than  +  40  nT  for  our  test 
subset.  Problems  with  magnetic  background  were  noted  by  both  the  ORNL  and  NRL  teams  that 
flew  surveys  of  this  area.  Again,  the  data  proved  to  be  multiffactal,  although  in  this  case  an 
isotropic  simulation  sufficed. 

Limitations  of  this  Pilot  Study 

One  major  challenge  we  faced  in  multifractal  characterization  of  the  airborne  magnetic  datasets 
collected  at  all  three  sites  was  determining  the  spectrum  and  multifractal  parameters  for  the 
magnetic  field  produced  by  the  geologic  background  in  the  presence  of  UXO-related  magnetic 
anomalies.  At  each  site  we  selected  portions  of  the  full  data  set  with  low  concentrations  of  UXO, 
and  in  some  cases  fitted  and  removed  UXO  anomalies,  but  the  possibility  remains  that  the 
magnetic  signatures  of  small  UXO  items  were  incorporated  in  our  characterization  of  geologic 
noise.  Ideally,  we  would  like  to  see  future  investigations  of  UXO  sites  include  collecting 
magnetic  data  over  nearby  regions  with  similar  geology,  but  no  history  of  ordnance  activity.  This 
would  improve  characterization  of  the  geologic  noise,  but  obviously  would  also  increase  the  cost 
of  data  acquisition.  However,  the  mobilization  and  demobilization  costs  would  not  change,  so 
the  incremental  cost  of  collecting  background  data  would  be  modest  and,  we  believe,  well  worth 
the  effort,  particularly  for  sites  with  potentially  high  geologic  noise  where  data  quality  assurance 
is  especially  import. 

Secondly,  because  most  contractors  have  focused  their  data  processing  stream  on  producing  an 
optimal  UXO  anomaly  map,  the  data  we  received  were  not  processed  for  total  magnetic  field. 
Both  ORNL  (later  this  team  moved  to  Battelle)  and  NRL  used  low-pass  filtering  to  reduce 
instrument  noise,  helicopter  noise,  geologic  background  and  other  noise  sources.  Instead  of 
producing  total  field  maps,  these  researchers  mapped  residual  magnetic  field  (e.g..  Figure  14), 
analytic  signal,  or  lists  of  anomaly  targets.  Because  low  pass  filtering  does  not  remove  geologic 
anomalies  with  wavelengths  comparable  to  the  UXO,  we  think  a  better  approach  is  to  understand 
the  full  spectrum  of  the  geologic  noise,  and  for  this  we  need  the  true  total  field  data.  As  a 
consequence,  we  had  to  reprocess  the  raw  data  provided  by  the  contractors  to  eliminate  heading 
error,  rotor  noise,  etc.,  to  produce  total  field  data.  Unquestionably  the  contractors,  who  are  more 
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Synthetic  Noise  with  25  Embedded  Targets 


Easting  -  (m) 

Figure  29:  Simulated  magnetic  background  (nanoTeslas)  and  dipole  UXO  targets.  The  actual 
locations  of  the  targets  are  marked  with  blue  circles.  The  anomalies  identified  by  thresholding 
are  marked  with  black  asterisks.  Numerous  false  positives  and  false  negatives  are  evident. 


cognizant  of  all  the  field  data  acquisition  and  processing  steps,  could  have  produced  better  total 
field  data  had  they  been  tasked  to  do  so  as  part  of  their  original  surveys. 

A  third  limitation  is  the  simplifying  assumption  that  magnetic  background  could  be 
approximated  as  coming  entirely  from  thin  surface  layer,  neglecting  the  distribution  of  magnetic 
susceptibilities  as  a  function  of  depth.  Pragmatically,  this  approximation  can  be  justified  because 
the  magnetic  field  above  the  surface  of  the  earth  created  by  any  3-D  distribution  of  magnetic 
susceptibilities  can  be  treated  as  coming  from  an  equivalent  source  layer  at  the  surface  (Blakely, 
1996).  Using  a  multifractal  representation  of  an  equivalent  source  layer  in  no  way  limits  the 
simulations  of  the  magnetic  field.  But  because  the  equivalent  source  layer  may  bear  little 
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resemblance  to  the  true  3-D  distribution  of  magnetic  susceptibilities,  it  would  be  more  rigorous 
to  use  a  full  3-D  multiffactal  representation  of  magnetic  susceptibility,  an  endeavor  that  we  wish 
to  pursue. 

Finally,  as  discussed,  our  current  methodology  assumes  isotropy,  or  an  anisotropy  that  can  be 
modeled  by  simple  stratification.  The  data  set  from  Fort  Ord  demonstrated  that  isotropy  is  not 
always  a  valid  assumption,  and  we  were  fortunate  that  in  this  case  a  simple  representation  of 
anisotropy  produced  simulations  matching  the  character  of  the  data  at  this  site.  However, 
developing  an  implementing  a  more  general  approach  will  be  a  priority  in  our  future  work. 

Potential  Applications 

This  pilot  study  focused  on  the  developing  methods  for  multiffactal  characterization  and 
simulation  of  airborne  magnetic  data.  The  next  steps  include  working  on  overcoming  the 
limitations  discussed  above,  and  on  application  of  the  simulations  to  quality  assurance  (QA)  for 
UXO  detection  and  discrimination.  To  give  an  indication  of  how  this  work  might  progress,  we 
created  a  simple  example.  Using  the  multiffactal  parameters  we  obtained  for  the  SI  target  area  at 
Pueblo  of  Isleta,  we  simulated  a  0.5  x  0.5  magnetic  field  at  a  flight  height  of  2.0  m,  with  25  UXO 
targets  embedded  at  random  x-y  locations,  buried  at  depths  of  0  to  1.5  m.  The  targets  were 
approximated  as  simple  dipoles  with  purely  induced  magnetic  moments,  and  varied  in  strength 
from  peak  amplitudes  of  8  -  150  nT.  Consequently,  some  of  the  targets  produced  anomalies 
clearly  above  the  magnetic  background,  others  were  buried  in  the  noise.  We  then  imported  the 
this  dataset  into  UXOlab  and  used  thresholding  at  35  nT  to  automatically  locate  anomalies. 

Figure  29  compares  the  actual  target  locations  with  the  target  picks.  At  this  threshold,  UXOlab 
pick  23  targets.  Nine  are  correctly  identified,  12  are  false  positives,  and  there  are  14  false 
negatives.  Clearly,  statistics  can  be  improved  by  manually  eliminating  picks  clustering  on  the 
same  peak,  testing  different  thresholds,  and  using  other  peak  picking  algorithms.  Furthermore, 
more  sophisticated  models  of  UXO  items  can  be  using  in  the  forward  simulation.  The  point, 
however,  is  that  a  realistic  model  for  the  geologic  background  facilitates  “what-if  ’  scenario 
testing  and  any  scale,  and  the  best  choice  of  algorithms  for  a  given  site. 

Future  Work 

The  ability  to  create  of  realistic,  site-specific  data  simulations  allows  the  fine-tuning  of  detection 
and  discrimination  algorithms.  The  key  next  step  is  to  develop  collaborations  with  SERDP 
researchers  working  on  these  algorithms  and  customize  our  simulations  to  their  needs.  Some 
customizations  would  be  straightforward,  such  as  simulating  measurements  collected  along 
irregular  flight  paths  and  at  uneven  altitudes,  or  simulating  vertical  magnetic  gradient  data. 

Others  enhancements  would  require  further  research. 

We  have  already  discussed  the  need  to  incorporate  a  generalized  anisotropy  in  the  multiffactal 
simulations.  To  do  this  also  requires  improved  data  characterization,  going  beyond  the  double 
trace  moment  method  and  isotropic  fitting  of  the  structure  function.  A  number  of  approaches  to 
extracting  anisotropic  multifractal  parameters  have  been  discussed  in  the  literature  (e.g.,  Lewis, 
et  al.,  1999;  Tennekoon  et  al.,  2003;  Cheng,  2004),  and  these  need  to  be  investigated  for  this 
application. 
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UBC  researchers  have  found  that  one  of  the  more  reliable  methods  of  distinguishing  between 
natural  and  UXO-related  anomalies  is  to  examine  the  orientation  of  the  anomaly's  dipole  field 
(e.g.,  Billings,  2004).  UXO  are  more  likely  to  have  dipole  moment  oriented  differently  than  the 
inducing  field.  To  test  this  we  would  need  to  refine  our  magnetic  model  to  allow  for  the 
inclusion  of  remnant  magnetism  in  the  simulated  geologic  background. 

Finally,  we  would  like  to  point  out  that  the  multifractal  simulation  approach  we  have  applied  in 
this  pilot  study  is  not  limited  to  the  simulation  of  magnetic  susceptibility.  It  is  feasible  to 
develop  simulations  for  physical  properties  that  govern  other  commonly  used  UXO  detection 
methods  such  as  electrical  conductivity  (electromagnetic  methods)  or  dielectric  properties 
(ground  penetrating  radar).  The  study  of  correlated  multifractal  fields  is  a  relatively  new 
research  area.  Marsan  et  al.  (1996)  discussed  the  theory  for  correlated  multifractal  fields,  but  we 
are  not  aware  of  any  work  that  addressed  the  problem  in  an  applied  sense.  Development  of  joint 
multifractal  simulation  of  multiple,  linked  geophysical  quantities  would  facilitate  the  current 
trend  toward  data  fusion  approaches  to  UXO  detection  and  discrimination. 
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